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ABSTRACT 


The  current  status  of  computational  capabilities  for  calculating  viscous- 
in viscid  transonic  Hows  other  than  the  solution  of  Navier-Stokes  equations  is 
presented.  Techniques  for  solving  transonic  inviseid  flows  and  compressible 
integral  boundary  layer  methods  are  reviewed,  and  systems  for  strong  viscous- 
inviscid  interactions  are  described.  Generally,  the  transonic  viscous-inviscid 
interaction  is  characterized  by  a  subcritical  boundary  layer  with  a  supersonic 
outer  stream.  The  thickening  boundary  layer  produces  a  pressure  rise  which 
causes  further  thickening  of  the  boundary  layer.  The  physical  flow  is  best 
modeled  by  direct  coupling  of  the  viscous  and  inviseid  systems  to  allow 
immediate  interaction  between  the  shock  wave  and  the  boundary  layer.  It 
appears  that  the  method  of  integral  relations  for  the  outer  inviseid  flow, 
combined  with  an  integral  boundary  layer  scheme,  possesses  such  a  capabil¬ 
ity.  To  facilitate  the  computation,  an  hybrid  approach  to  the  transonic 
inviseid  solution,  which  consists  of  the  finite  difference  method  for  solving 
the  overall  transonic  inviseid  potential  flow  field  and  the  method  of  integral 
relations  for  solving  Euler’s  equation  in  the  shock  region,  is  suggested. 
Finally,  the  application  of  the  steady  two-dimensional  methods  to  the  quasi 
two-dimensional  problem  on  axisymmetric  stream  surface  of  a  cascade  flow 
at  transonic  speeds  is  discussed. 
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I.  INTRODUCTION 


Recent  advances  in  numerical  techniques  have  provided  a  powerful  tool  in  transonic  How 
research.  A  complete  full  potential  How  solution  is  now  attainable  in  a  matter  of  minutes  by  the 
finite  difference  relaxation  scheme,  liven  a  fully  inviscid  (nonisentropic)  solution  can  be  deter¬ 
mined  numerically  by  the  unsteady  finite  difference  approach  at  the  cost  of  large  computer  time 
or  by  the  method  of  integral  relations  with  certain  iterative  procedures.  The  real  transonic  flow 
problem,  however,  is  complicated  by  the  viscous  effect.  The  strong  viscous-in  viscid  interaction 
caused  by  the  shock  wave  thickens  the  boundary  layer  rapidly,  and  the  flow  eventually  separates 
from  the  surface.1'6  In  such  cases,  the  surface  flow  properties  cannot  be  specified  in  advance,  as 
in  the  usual  formulation  of  boundary-layer  theory,  but  are  determined  by  interaction  of  the 
outer  inviscid  flow  and  the  inner  viscous  flow  near  the  surface.  This  has  been  referred  to  as  the 
transonic  viscous-in  viscid  interaction  problem. 

The  problem  of  the  viscous-inviscid  interaction  in  the  transonic  region  is  complicated  by  the 
fact  that  the  inviscid  flow  field  is  governed  by  mixed  elliptic  and  hyperbolic  partial  differential 
equations  for  a  compressible  fluid  flow.  Because  the  major  portion  of  the  flow  field  is  of  the 
elliptic  type,  the  velocity  at  the  edge  of  the  boundary  layer  at  any  location  depends  on  the  com¬ 
plete  displacement  thickness  distribution.  At  the  same  time,  the  solution  must  satisfy  the  con¬ 
straints  characterized  by  the  inviscid  supercritical  flow  the  regularity  condition  at  sonic  points, 
fhe  complete  solution,  therefore,  involves  tedious  iterative  procedures.  Some  advances  in  the 
area  of  laminar  viscous-inviscid  interactions  at  transonic  speeds  were  made  by  MacCormack,7 
Klineberg  and  Steger,8  Brilliant  and  Adamson,9  and  Tai.10  Although  the  laminar  flow  yields  an 
adequate  model  for  assessing  the  basic  mechanism  of  the  strong  viscous-inviscid  interaction  proc¬ 
ess,  Hows  of  practical  interest  at  transonic  speeds  are  turbulent  because  of  high  Reynolds  number 
conditions.  Theoretical  treatment  of  the  transonic  turbulent  viscous-inviscid  interaction  is  fur¬ 
ther  complicated  by  the  turbulence.  In  the  past  several  years,  remarkable  efforts  have  been 

1  1  IQ 

devoted  to  the  topic.  The  problem  has  been  treated  with  various  levels  of  assumptions.  In 

most  cases,  however,  theoretical  considerations  were  limited  to  weak  shock  conditions.  The 
problem  of  strong  viscous-inviscid  interaction  at  transonic  speeds  was  discussed  by  Klineberg  and 
Steger  8  and  Tai. 10,15 

Generally,  the  problem  of  transonic  viscous-inviscid  interaction  is  best  described  by  the 
Navier-Stokes  equations.  Because  of  extensive  computer  capacity  requirements  and  other  nu¬ 
merical  problems,  a  complete  Navier-Stokes  solution  for  the  subject  problem  still  seems  to  be 
some  time  in  the  future.19  Researchers  will  continue  for  many  years  to  rely  on  other  procedures 
to  adequately  account  for  the  viscous  effects.  Among  various  techniques,  the  coupling  of  a  valid 
transonic  inviscid  solution  with  an  integral  boundary-layer  method  seems  to  offer  an  adequate 
approach  to  the  problem. 

In  the  present  paper,  therefore,  a  review  of  both  transonic  inviscid  solutions  and  the  integral 
boundary-layer  methods  will  be  presented  in  Sections  II  and  IN.  Because  the  review  will  not  be 
exhaustive  to  cover  such  a  field  of  phenomenal  research  activities,  emphasis  will  be  placed  on 
prediction  methods,  either  viscous  or  inviscid,  which  are  adequate  for  handling  the  subject  of 


viscous-in  viscid  transonic  Hows.  In  Section  IV,  the  problem  of  transonic  viscous-inviscid  inter¬ 
action  is  described  and  calculation  procedures  are  presented  along  with  some  theoretical  results. 
Finally,  the  application  of  the  methods  developed  for  Hows  around  airfoil  to  cascade  How  prob¬ 
lems  is  discussed. 


II.  REVIEW  OF  TRANSONIC  INVISCID  FLOW  SOLUTIONS 

Theoretical  methods  for  solving  the  transonic  inviscid  How  have  been  advanced  substantially 
in  the  past  ten  years.  Because  of  the  nonlinear  nature  of  the  tlow  equations,  virtually  all  the 
solution  procedures  have  been  developed  numerically.  Solutions  are  achievable  only  with  the  aid 
of  high  speed  computers. 

1.  BASIC  FLOW  EQUATIONS 

The  basic  system  of  equations  that  govern  an  inviscid  compressible  fluid  flow,  without  body 


force,  is: 

Continuity 
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where 

p  =  static  density,  P  =  static  pressure,  q  =  velocity, 
t  =  time,  S  =  entropy,  T  -  temperature, 

e  =  internal  energy. 

Eq.  (3)  can  be  used  to  compute  the  rate  of  change  of  entropy  of  a  fluid  particle  associated  witli 
the  shock  wave  in  a  supercritical  flow. 

For  regions  prior  to  the  shock  wave  and  the  far  Held,  the  changes  of  state  of  the  fluid  par¬ 
ticle  are  isentropic 


DS, 

TT"~  =  0  or  Sj  =  Constant 


(4) 
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Instead  of  the  energy  equation,  therefore,  a  simpler  isentropic  relation  may  be  used 


(5) 


For  regions  following  the  shock  wave,  where  the  fluid  has  a  finite  increase  in  entropy  (or  a 
decrease  in  the  total  pressure),  the  isentropic  relation  applies  along  a  streamline 
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where  subscripts  1  and  2  denote  states  before  and  after  the  shock  wave,  and  o,  the  stagnation 
value.  ev  is  specific  heat  at  constant  volume. 

The  change  of  entropy  (S->- S|  )  is  normally  small  for  transonic  Hows.  However,  because  of 
the  corresponding  change  of  the  downstream  conditions,  it  may  influence  the  location  ol  the 
shock  wave  and  thus  alter  considerable  portions  of  the  entire  flow.  Nonetheless,  the  flow  behind 
the  shock  wav'e  is  still  isentropic  along  a  streamline.  The  new  isentropic  relation,  however,  must 
be  based  on  a  new  entropy  level  which  is  slightly  higher  than  its  freestream  value.  The  new 
entropy  level  differs  from  one  streamline  to  another  because  the  shock  strength  encountered  on 
each  streamline  differs.  Therefore, 


S',  ¥=  Constant 


(7) 


Because  of  the  current  status  of  numerical  solutions  for  transonic  viscous-in  viscid  inter¬ 
action  Hows,  the  present  paper  will  cover  only  two-dimensional,  steady  flow  conditions.  A  brief 
discussion  on  the  application  to  the  unsteady,  three-dimensional  flow  associated  with  turbo¬ 
machinery  will  be  given  in  Section  V. 

a.  Full  Inviscid  Flow  Equations 

For  two-dimensional,  steady  Hows,  the  preceding  system  of  equations  reduce  to  the  follow¬ 
ing  form  in  Cartesian  coordinates  (x,  y): 
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where  U  and  V  are  velocity  components  and  note  that  S|  =  constant  and  S->  =£  constant. 


b.  Full  Potential  Flow  Equation 


It*  the  entire  How  field  is  assumed  isentropic,  the  flow  is  everywhere  irrotational.  The  veloc¬ 
ity  is  related  to  a  potential  function.  <1>.  as 


d<I> 
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(  12) 


Phe  exact  equation  for  the  potential  function  for  two-dimensional  compressible  How  can 
be  written  in  Cartesian  coordinates  as 


where 
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and  a  is  the  speed  of  sound.  Eq.  ( 13)  is  elliptic  if  the  local  Mach  number  M  <  I  and  hyperbolic 
if  M  >  1 .  t  hus,  Eq.  (8)  through  (11)  are  reduced  to  a  single  differential  equation.  Eq.  ( 13). 
Once  the  solution  is  found,  the  velocities  are  calculated  from  Eq.  ( 12).  Eq.  ( 13)  is  valid  only  in 
cases  with  weak  shock  waves  where  the  change  of  entropy  can  be  neglected.  By  definii  <v.  of  1 
irrotational  flow,  there  is  no  transonic  wave  drag. 

c.  Small  Disturbance  Equation 

Further  simplification  of  the  equation  system  is  to  assume  the  perturbation  velocities  being 
small,  resulting  in  the  transonic  small  disturbance  equation 
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where  0  is  the  perturbation  potential  function  and  other  symbols  have  usual  meaning. 

The  validity  of  t his  equation,  of  course,  bears  directly  on  the  assumptions  made,  i.e.,  weak 
shock  waves,  thin  airfoil  at  small  angle  of  attack.  And  again,  there  is  no  transonic  wave  drag. 

In  all  the  foregoing  systems,  the  viscosity  is  commonly  neglected.  The  problem  in  coupling 
the  viscous  solution  will  be  discussed  later. 

2.  NUMERICAL  SOLUTIONS 

I  or  convenience  of  discussion,  numerical  methods  for  the  transonic  inviscid  problem  can 
be  classified  into  (a)  finite  difference  method,  (b)  method  of  integral  relations,  (c)  hybrid  meth¬ 
od.  id)  finite  element  method,  and  (e)  integral  equation  method.  A  brief  review  of  these  meth¬ 
ods  i>  presented  below. 

a.  Finite  Difference  Method 

The  finite  difference  techniques  tor  solving  the  transonic  flow  equation  were  introduced 
earlier  by  Emmons  in  lc)48“°  using  central  differences  in  the  subsonic  region  and  upwind  differ¬ 
encing  in  the  supersonic  /.one.  The  technique  was  greatly  improved  by  Murman  and  Cole-1  along 
with  a  successive  line  over-relaxation  scheme.  Through  subsequent  improvements  in  the  area  of 
different  operators—  and  suitable  and  faster  relaxation  schemes the  method  has  been 
highly  developed  for  solving  the  full  potential  flow  equations,  as  well  as  the  transonic  small- 
disturbance  equations.  Survey  papers  given  by  Yoshihara-9  and  Bailey30  summarizing  early 
work  in  the  use  of  the  finite  difference  method  and  those  of  Jameson,31  Murman,3-  and 
Ballhaus  33  cover  techniques  implemented  in  the  current  widely  used  codes  such  as  those  by 
Bauer  et  al./4  Jameson,-5  and  Carlson.35 

To  illustrate  the  finite  difference  scheme  in  solving  the  transonic  potential  flow  equation,  a 
comprehensive  discussion  is  given  by  Murman.3-  For  simplicity,  however,  the  approach  is  easily 
understood  with  reference  to  a  Cartesian  coordinate  system.  For  this  reason,  the  analysis  devel¬ 
oped  by  Carlson35  using  Cartesian  coordinates  is  summarized  below. 

The  full  potential  flow  equation,  F.q.  (13),  can  be  recast  into  the  form  of  a  perturbation 
potential  0 


(a-  -  U- )  0XX  -  2UV0xy  +  (a:~V-)0yy  =  0 


(16) 


where  U  and  V  are  given  by 

U  =  =  q^  (cos  a  +  0X ) 

v  =  =  Moo  (sin  a  +  0y  > 
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along  with  boundary  conditions  of  velocity  normal  to  the  airfoil  surface  being  zero  and  a  veloc¬ 
ity  potential  satisfying  field  behavior  at  infinity. 

To  incorporate  the  coordinate  stretching  in  the  physical  plane  and  to  avoid  computational 
problems  in  the  supersonic  region,  the  potential  How  equation  is  further  arranged  as 


~  a^)^SS  +  ^NN  ~  ° 

where 

^SS  =  ~2  [u2f<typ|  +  2UVfg*iT?  +  V2g(g<AT?)J?j 
0NN  =  72  [v2f(f*t>{  -  2UVf^n  +  U2g(g0TJ)J 


(18) 


(19) 


The  S  and  N  denote  the  vector  parallel  and  perpendicular  to  the  local  velocity,  respectively,  and 
f  and  g  are  functions  relating  the  Cartesian  coordinates  (x,  y)  and  stretched  coordinates  (£.17), 
respectively,  as 


d£  dr? 

dx  '  8  ~  dy 


CO) 


The  finite  difference  expressions  with  reference  to  a  mesh  system  shown  in  Fig.  1  are,  for 
example,  for  contributions  to  0^  when  q-  >  a~  : 


(,<V?  ^{‘i  +  'A  (^i  +  1  ,  j  ~  *i  -  !*  (^ij  1 ,  j*| 

For  contributions  to  0s<j  when  q-  >  a^  and  V  >  0 
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These  equations  are  so-called  upwind  difference  formulas  for  supersonic  points.  For  the 
subsonic  region,  a  central  differencing  scheme  is  employed,  for  example. 


where  the  relaxation  factor  w  lias  been  incorporated  into  the  difference  formulas.  When  these 
expressions.  Kqs.  (21)  through  (23),  arc  substituted  appropriately  into  Eqs.  ( 18)  and  ( 19)  or 
Fq.  { lb),  the  result  is  a  tridiagonal  system  of  equations  that  can  be  solved  for  the  current  values 
of  the  function  0  column  i.  Numerical  stability  is  achieved  by  the  incorporation  of  both  old  0 
and  new  0+  values  into  the  finite  difference  formula  and  the  inclusion  of  time-like  derivatives 
in  the  relaxation  process.  The  result  obtained  by  Carlson’s  code77  is  very  close  to  that  using 
the  Jameson  coder''  see  Fig.  2.  Other  applications  of  the  finite  difference  method  were  re¬ 
ported  elsewhere.7*1  ,77 

An  important  feature  in  formulation  of  difference  equations  is  the  consistency  condition 
with  the  original  partial  differential  equations.  Murman77  discussed  the  topic  in  great  detail  and 
classified  the  difference  equations  for  transonic  flows  into  fully  conservative  form  and  not  fully 
conservative  form  (or  quasilinear  form).  Fig.  3  shows  the  difference  in  computed  results  using 
these  two  forms.  It  is  observed  that  the  nonconservative  form  yields  a  weaker  shock  wave  but 
located  slightly  more  upstream  than  that  determined  by  the  fully  conservative  equations. 

Murman  has  pointed  out  deficiencies  in  using  the  nonconservative  formulas.  Nonetheless,  the 
nonconservative  system  is  preferred  by  most  users  because  of  Its  better  correlation  with  the 
experimental  data.  Tai  remarked  that  the  momentum  deficiency  resulting  from  the  non¬ 
conservative  form  has  similar  effects  on  the  strength  and  location  of  the  shock  wave  resulting 
from  the  total  pressure  deficiency  aft  the  shock  caused  by  the  entropy  rise  across  the  shock  wave. 
It  will  be  further  discussed  in  the  next  section. 

b.  Method  of  Integral  Relations 

The  method  of  integral  relations  (MIR)  provides  another  avenue  to  the  solution  of  transonic 
How  problems.  Previous  applications  of  the  method  to  transonic  flows  past  airfoils  include  those 
by  Holt  and  Masson,79  Melnik  and  Ives,40  Sato,41  and  Tai4^  for  various  How  conditions  in  the 
transonic  regime.  A  review  of  the  method  in  transonic  flow  is  given  by  Holt  47  The  main  advan¬ 
tage  in  using  MIR  is  its  ability  to  solve  the  full  inviscid  How  equations  directly;  thus  the  assump¬ 
tion  of  isentropic  flow  is  not  necessary. 

In  the  numerical  solution,  it  is  convenient  to  write  equations  in  nondimensional  form 
normalized  by  freestream  values.  Therefore,  the  full  inviscid  How  equations,  (F.qs.  (8)  through 
(II)),  are  recast  as 
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Continuity 

i 

x-Momentum 


y-Momentum 

i 


where 


3  (pU)  3  (pV ) 
3x  +  3y 


(24) 


—  (KP  +  pU:)  +  2.<pUV)  =  0 
dx  dy 


^(pUV)  +  —  (KP  +  pV- )  =  0 


K  =  1/(7M£). 


(25) 


(26) 


(27) 


The  boundary  conditions  are  as  follows:  at  the  airfoil  surface,  the  normal  velocity  com- 
i  ponent  equals  to  zero,  i.e.. 


! 


J 

J 


% 


=  0 


and  at  infinity,  the  flow  is  undisturbed,  i.e,, 

P  =  p  =  U  *  I 
V  =  0 

Briefly,  in  applying  the  system  of  flow  equations  must  be  written  in  divergence  form 

3  a 

^A(x,y,U - )  +  ^B(x,  y,  U - )  =  Q(x,y,U)  (28) 

The  divergence  form  of  Hqs.  (24)  through  (26)  may  then  be  integrated  outward  from  the  airfoil 
surface  (but  not  necessarily  normal  to  the  surface)  to  each  strip  boundary  successively  at  some  x 
station.  This  procedure  reduces  the  partial  differential  equations  to  ordinary  ones.  To  perform 
the  integration,  the  integrand  is  approximated  by  interpolation  polynomials,  for  example  A,  by 


A 


N 

£  ak  <x)(y  -  y0)k 
k  =  0 


(29) 
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v\  here 


R*  =  the  number  of  strips 
a^t  x)  =  constants  evaluated  at  strip  boundaries 
\0  =  the  location  of  the  base  strip  boundary. 

In  principle,  the  actual  flow  variation  may  be  represented  more  closely  by  an  increasing  number 
of  strips. 


Using  a  basic  second-order  approximation  for  Hq.  <2C>),  the  method  can  be  implemented 
with  up  to  five  or  six  strips  for  desired  accuracy.  The  process  is  illustrated  in  Fig.  4.  The  idea  is 
to  treat  the  whole  integration  domain  as  a  series  of  different  effective  regions;  a  small  number  of 
strips  may  be  used  within  each  region.  Three  effective  regions  are  designated  by  strip  boundaries 
(0.  1 .  2f  (a,  b.  c>.  and  (£,  77.  ?).  In  each  effective  region,  the  How  tie  Id  is  divided  into  two  strips 
and  approximated  in  the  usual  way  by  a  second-order  polynomial.  The  MIR  is  first  applied  to 
the  boundaries  £,  77.  f  with  the  purpose  of  determining  the  shape  of  a  streamline  r\  and  tlow  con¬ 
ditions  along  it.  On  the  uppermost  boundary  f.  free  stream  conditions  are  assumed  to  apply. 
Secondly,  MIR  is  applied  along  boundaries  a,  b,  and  c  to  determine  conditions  along  a  streamline 
b,  closer  to  the  profile.  Finally,  MIR  is  applied  in  the  disturbed  part  of  the  field  along  the 
boundaries  0,  1 ,  and  2. 

The  resulting  ordinary  differential  equations  take  on  the  form  along  boundaries  a  and  b.  tor 
example 


—  =  f,  <x.y.Ua.Va....) 


—  =  t:,x.y.Ua.Va....> 


dUh 

dx 


g,(x.y.Ua.va. ...» 


l'vh 

—  =  p:(x.y.Ua.Va....) 


v  c’-1  / 


= 


( 30* 


(31 1 


(32) 


(33) 


(34» 
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(35) 


•e 


w  Ik*  iv 


<  3f> ) 


(37> 


I  Ik*  set  ot  Fqs.  (31 )  through  (34)  are  then  solved  hy  the  fourth-order  Runge-kutta  scheme. 
I  he  solution  is  obtained  hy  initial  value  techniques  associated  with  a  two-point  boundary-value 
problem  For  t he  flow  over  a  lifting  airfoil,  the  complete  solution  procedure  consists  of  iterative 
processes  lor  handling  the  regularity  condition  at  the  saddle  point  (sonic  point),  for  determining 
the  shock  location  and  for  enforcing  the  Kutta  condition. 

I  he  method  allows  the  exact  Rankine-Hugoniot  relations  to  be  applied  at  the  shock  tor 
determining  the  flow  properties  aft  of  the  shock 

( 7  +  1  )  M  | 

p:  =  pi  (7  -  i  >  M ,  +2 


Pi 


f>l 


2y  <M“  -  1  ) 

7  +  1 


The  change  of  entropy  is  obtained  in  terms  of  Pi  and  p-,: 


Si 


( 3*> ) 


(40) 


where  Si  varies  from  one  streamline  to  the  other. 

Figs.  5  and  (»  show  typical  results  obtained  by  the  method  of  integral  relations.  These 
figures  are  taken  from  Tai4*  and  include  data  from  Magnus  and  Yoshihara  44  Stivers  4> 
(iarabedian  and  Korn  40  Krupp  and  Murrnan  47  and  Kacpriynski  et  al  48  Reasonable  agreements 
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between  the  theory  and  experiment  were  achieved  except  in  the  neighborhood  of  the  shock 
wave.  The  disagreement  there  is  attributed  to  a  viscous-in  viscid  interaction  between  the  shock 
wave  and  boundary  layer.  Note  the  effects  of  change  of  entropy  (Fig.  5)  and  of  the  shock  angle 
(t  ig.  (>)  were  explicitly  evaluated.  The  increase  in  entropy  corresponds  to  the  decrease  in  total 
pressure  4t)  Physically,  the  change  of  entropy  creates  vorticity  behind  the  shock  wave;  it  has  a 
cumulative  effect  in  the  far  downstream  and  consequently  feeds  back  to  the  shock  itself  4~  In 
the  case  of  a  finite  increase  in  entropy  across  the  shock  wave,  it  was  found  that  the  shock 
strength  was  weakened  and  the  shock  location  moved  forward.  It  has  a  similar  effect  of  the 
momentum  deficiency  of  the  nonconservative  form  in  the  finite  difference  method,  although 
their  mathematical  bases  are  completely  different. 

c.  Hybrid  Method 

The  finite  difference  scheme  and  the  method  of  integral  relations  for  solving  the  transonic 
flow  problems  offer  different  kinds  of  advantages  and  disadvantages.  The  former  is  well  devel¬ 
oped,  and  easy  for  the  user  to  implement;  it  is  limited  to  the  isentropic  How  assumption  (weak 
shock  wave  assumption)*  and  suffers  mathematical  inconsistency  of  the  nonconservative  formu¬ 
la.  For  the  method  of  integral  relations,  on  the  contrary,  the  situation  is  opposite.  It  solves 
Filler's  equation  with  a  small  computation  requirement;  but  because  of  the  multiple  iterative 
processes  involved,  the  whole  solution  procedure  cannot  be  automated  without  man-machine 
interactions.  To  take  advantage  of  the  both  and  yet  avoid  their  shortcomings,  a  hybrid  method 
combination  of  the  finite  difference  method  and  the  method  of  integral  relation  appears  to 
offer  an  unique  capability  for  solving  transonic  flow  problems  with  shock  waves. 

The  central  idea  for  the  hybrid  method  proposed  by  Tui,  is  illustrated  in  Fig.  7.  The  overall 
mixed  fiowfield  is  governed  by  the  potential  tlow  equation  except  the  shock  region.  In  the  shock 
region  and  the  region  downstream  of  the  shock,  the  How  is  governed  by  Euler's  equation  that 
allows  entropy  change,  creating  the  transonic  wave  drag.  First,  the  potential  tlow  equation  is 
solved  by  the  finite  difference  relaxation  procedure.  The  flow  properties  so  calculated  are  valid 
up  to  the  boundaries  of  the  shock  region  marked  by  solid  dots  in  Fig.  7.  The  values  at  these 
boundaries  are  then  taken  as  the  initial  condition  (those  along  the  normal  line)  and  the  boundary 
condition  (those  along  the  streamwise  line)  for  the  numerical  solution  by  the  method  of  integral 
relations.  In  so  doing,  the  major  difficulty  in  treating  the  saddle  point  (sonic  point)  regularity 
condition  associated  with  the  method  of  integral  relations  is  avoided  and,  therefore,  in  general, 
the  procedure  can  be  fully  automated. 

The  concept  is  now  being  developed  into  a  complete  numerical  procedure. 


♦The  finite  difference  approach  using  Euler’s  equations  has  been  developed  using  unsteady  techniques44  with 
extensive  computer  time  and  storage  requirements. 
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d.  Finite  Element  Method 

The  use  of  the  finite  element  method  in  transonic  Hows  is  relatively  new  A  lew  references 
can  be  found  in  the  literature. 1  Because  of  problems  of  hyperbolic  behavior,  the  work  has 
been  limited  to  the  small  disturbance  equations. 

The  following  brief  discussion  of  the  method  is  taken  from  C  han  and  Brashears.S° 

In  using  the  finite  element  method,  the  transonic  small  disturbance  equation,  Eq.  ( 15).  is 
rewritten  as 

0Xx  +  ^yy  ~  ^ 

where 


f  -  (I+7)  0xj 


along  with  appropriate  boundary  conditions 


V0  =  o 


on  the  airfoil 
at  infinity 


Using  Galerkin's  weak  form,  Eq.  (41 )  is  then  transformed  to 


Jf* 


i  (Nj,xx  +  Nj,yy> tlxdy 


♦i  /ft 


fdxdy  =  0 


In  Eq.  (43):  Nj  and  Nj  represent  the  shape  functions.  0j's  are  the  unknown  parameters,  the 
indices  iand  j  run  from  I  to  the  total  number  of  unknown  parameters,  and  the  integrations  are 
to  be  performed  over  the  entire  domain  under  consideration. 

Upon  integrating  by  parts,  there  results 


Sjj0j  -  L, 


where 


sij  ~ffm l.«  N,.x  +  Ni.y 


Li  _  4  ■  /n  * ?• Ni  < Cl 


3 


1  lu‘  original  transonic  problem  is  tluis  transformed  into  a  system  of  algebraic  equations  in  the 
form 

Sjj  *j"‘  =  L‘n-'>  (47) 

w  here  0,n  1  is  the  n1*1  step  solution  and  L(n  *  1  *  is  the  load  matrix  evaluated  by  using  the  previous 
solution  l*ii .  (47)  is  to  be  solved  subject  to  certain  prescribed  convergence  criterion. 

I  ie  S  shows  a  sample  result  of  the  finite  element  method  for  a  supercritical  flow  past  an 
VAC  A  MAOOo  airfoil;  taken  from  Kef.  50.  Using  a  triangular  cubic  element  with  136  modes, 
the  agreement  between  theory  and  experiment^-  is  only  qualitative.  Generally,  the  method 
umlerpreiliets  the  flow  before  the  shock  wave.  Comparison  between  the  finite  element  solution 
and  the  finite  difference  solution  of  a  transonic  flow  is  made  by  Hafez  et  al.,^*  as  shown  in  Fig.  9. 
Again,  the  finite  element  method  seems  to  underpredict  the  flow  near  the  leading  edge  of  the  air- 
tod  Various  ways  for  improving  the  finite  element  method  in  handing  the  mixed  flow  region 
and  extension  to  the  full  potential  llow  equations  have  been  discussed  in  Rets.  50  and  51 . 

e.  Integral  Equation  Method 

As  opposed  to  the  aforementioned  methods,  the  integral  equation  method  has  been  re¬ 
garded  as  a  semmumerical  method  in  the  sense  that  it  employs  analytical  functions  for  represent¬ 
ing  the  integrand.  I  he  integral  equation  was  derived  by  Oswatitsch*^  thirty  years  ago  for  the 
solution  of  tin*  transonic  small-disturbance  equation.  The  method  was  further  explored  by 
( i nl Kt rand  ,>4  Spreiter  and  Alksne.vS  and  Norstrud.S6  The  integral  equation  method  requires 
least  numerical  work.  The  disadvantage  of  the  method  is  the  difficulty  in  achieving  an  adequate 
field  integral  for  approximating  the  complicated  mixed  How  field.  Nixoir  •*  suggested  that  the 
accuracy  of  the  integral  equation  method  can  be  improved  considerably  by  applying  the  field 
integral  in  regions  of  strips  which  divide  the  entire  flow  field.  The  accuracy  of  the  integral  can 
be  increased  simply  by  increasing  the  number  of  strips.  The  computational  region  can  be,  there¬ 
fore .  extended  to  ‘infinity"  without  introducing  further  numerical  complexity.  Fig.  10  shows 
some  ty  pical  results  obtained  by  using  the  extended  integral  equation  method 
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III.  REVIEW  OF  INTEGRAL  BOUNDARY-LAYER  METHODS 


1.  LAMINAR  FLOWS 


a.  Laminar  Boundary-Layer  Equations 


The  governing  equations  for  a  compressible  laminar  boundary-layer  How  in  coordinates 
parallel  and  normal  to  the  surface  (s,  n)  are,  with  flow  properties  normalized  by  their  freest  ream 
values. 


Continuity 


3  (pu)  3  (pv) 
3s  +  3n 


(48) 


s-M  omen  turn 


3u  3u  dP  1  3  f  3ul 

3s  +  3n  K  ds  +  Poo Vqq c  3n  3n  J 


(49) 


where  c  is  the  chord  length  of  an  airfoil. 

The  boundary  conditions  are  as  follows:  at  the  surface,  u  =  v  -  0;  at  the  edge  of  the 
boundary  layer,  u  =  U  (s) 

b.  Klineberg-Lees  Method 

The  method  was  originally  developed  for  analyzing  the  supersonic  viscous-inviscid  inter¬ 
action  problem.59  Basically,  the  compressible  flow  equations,  Eqs.  (48)  and  (49)  are  trans¬ 
formed  into  an  equivalent  incompressible  form  by  means  of  the  Stewartson  transformation.60 
The  velocity  profile  based  on  an  incompressible  similarity  solution  can,  therefore,  be  used  in  the 
integral  approach.  The  resulting  partial  differential  equations  in  a  transformed  incompressible 
plane,  along  with  the  equation  for  the  rate  of  change  of  mechanical  energy,  i.e.,  the  moment  of 
momentum,  are  then  integrated  across  the  boundary  layer  from  n  =  0  to  5.  In  performing  the 
integration,  the  parameter  a  is  employed  to  denote  the  velocity  profile  for  both  attached  and 
separated  boundary  layers.  Profile  quantities  are  then  defined  as  functions  of  a. 

The  resulting  ordinary  differential  equations  are 


d5* 

—  =  F,  (8*.  a.  Mc)  <50> 


da 

d7 


F2  (5*.  a,  Me) 


(51) 
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i 

_ Li 


(52) 


dMe 

—  =  F3(S;,a.Me) 


where 


5*  =  boundary  layer  displacement  thickness 
a  =  profile  parameter 

M  =  Mach  number  at  the  edge  of  boundary  layer. 

I*  unctions  bp  Fi ,  and  depend  on  variables  6  •  ,  a,  Me  and  intermediate  parameters 
IF  J.  (J,  R.  and  Z,  which  have  been  defined  as  functions  of  a  in  Ref.  59.  Curve-fitted  poly¬ 
nomials,  based  on  the  similarity  solution  of  the  classical  boundary-layer  theory,  are  given  in  Ref. 
59  for  attached  and  separated  flows,  and  in  Ref.  61  for  wake-reverse  and  wake-forward  flows. 
The  original  ordinary  system  given  by  Ref.  59  has  a  temperature  parameter.  To  simplify  the 
discussion,  it  has  been  excluded  in  the  above  system  for  adiabatic  flows.10 

The  basic  features  of  the  method  are  that  ( 1 )  the  boundary-layer  properties  depend  on  the 
variable  a,  which  directly  describes  the  velocity  distribution,  and  (2)  the  velocity  at  the  edge  of 
the  boundary  layer  is  treated  as  a  dependent  variable,  rather  than  a  given  quantity  as  in  the  con¬ 
ventional  boundary-layer  approach.  The  static  pressure  can,  therefore,  be  determined  by  the 
interaction  between  the  outer  inviscid  flow  and  the  inner  viscous  layer  near  the  surface.  Such 
an  arrangement  allows  direct  coupling  of  the  viscous  system  with  inviscid  equations  for  handling 
the  viscous-inviscid  interaction  problem  to  be  discussed  later. 

Fig.  1  1  shows  the  results  of  a  typical  computation  for  an  adiabatic  flow  by  the  Klineberg- 
Lees  method.  The  two  solutions  represent  interactions  beginning  in  the  weak-interaction  zone 
for  both  a  finite  and  a  semi-infinite  flat  plate.  It  is  evident  that  a  small  change  in  the  perturba¬ 
tion  value  can  generate  very  different  integral  curves  because  of  the  high  degree  of  nonlinearity 
of  the  equations  in  the  subcritical  region.  The  weak  interaction  solution  corresponds  to  the 
conventional  boundary-layer  formulation,  i.e.,  the  Me  is  “impressed”  by  the  inviscid  solution 
rather  than  calculated  by  the  viscous-inviscid  interaction.  Note  that  the  boundary-layer  dis¬ 
placement  thickness  rises  drastically  due  to  the  strong  interaction. 

2.  TURBULENT  FLOWS 

a.  Turbulent  Boundary-Layer  Equations 

The  governing  equations  for  a  compressible,  adiabatic,  turbulent  boundary  layer  in  coor¬ 
dinates  parallel  and  normal  to  the  surface  (s,  n)  are  with  flow  properties  normalized  by  their 
tree  stream  values. 
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Continuity 
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where  j3  is  an  eddy  viscosity  parameter. 

The  boundary  conditions  are  as  follows:  at  the  surface,  u  =  v  =  0.  At  the  edge  of  the 
boundary  layer,  u  =  U  (s). 

b.  Kuhn-Nielsen  Method 

In  the  Kuhn-Nielsen  method  the  above  equations  are  transformed  into  incompressible 
plane  (£,  17)  with  the  aid  of  the  Stewartson  transformation60  along  with  the  assumptions  that  the 
viscosity  varies  linearly  with  the  temperature  and  that  the  Prandtl  number  equals  unity.  The 
velocity  profile  in  the  transformed  plane  used  in  performing  the  integration  across  the  boundary 
layer  is 


u  =  UT  |2.S  Cn(i  +7J+)  +  5.1  -  (3.397j+  +  5.l)e-°37,J+j 

-cosK)] 


(55) 
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where 


tl/3  ~ 


wake  velocity 
friction  velocity 
u  Tr\iv 


It  is  composed  of  an  inner  part,  consisting  of  a  laminar  sublayer  and  the  law-of-the-wall  function, 
and  an  outer  part,  a  wake  function.  It  is  a  modification  of  Cole's  law6*"  with  a  laminar  sublayer 
added. 

The  eddy  viscosity  model  used  in  the  Kuhn-Nielsen  method  is  an  extension  of  the  two- 
layer  model.  For  attached  How,  the  inner  layer  is  represented  by  an  expression  given  by 
Kleinstein.6^ 
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In  the  outer  layer,  the  ('la user  expression  along  with  an  interniittency  function  of  Klebanolf*4 
is  employed 
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For  separated  (low,  the  entire  profile  is  based  on  a  form  given  by  Alberft* 
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There  results  two  ordinary  differential  equations 
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where  An  through  A>3  are  defined  in  Ref.  12.  The  dependent  variables  are  ur.  6.  and  Ue.  In 
the  computation,  Ue  is  prescribed  except  in  the  region  near  separation,  uT  is  prescribed,  and  Ue 
is  calculated. 

The  results  of  a  typical  calculation  for  the  flow  over  a  bump  at  transonic  speed  are  shown  in 
Fig.  12.  The  viscous  and  inviscid  solutions  agree  within  one  percent.  For  comparison,  the  meas¬ 
ured  b*  and  separation  point  are  also  shown.  The  calculated  pressure  is  slightly  higher  than 
the  experimental  data.^  The  predicted  5*  is  in  excellent  agreement  with  the  data  upstream  ot 
the  shock  wave.  The  theoretical  b*  for  the  converged  strong  interaction  solution  downstream  of 
the  shock  is  slightly  lower  than  the  experimented  value  probably  due  to  the  uncertainty  in  spec¬ 
ifying  the  boundary-layer  conditions  downstream  of  the  shock.  More  thorough  analysis  ot  the 
local  shock -boundary  layer  interaction  is  needed  to  improve  this  aspect  of  the  calculation  and  to 
handle  the  case  of  separation  right  at  the  shock. 
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c.  Nash-Macdonald  Method 


The  Nash-Macdonald  method67  is  a  semi -empirical  formula  for  the  solution  of  the  momen¬ 
tum  thickness.  0.  in  the  integral  houndary-layer  equation  given  by  von  Karmun 
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where  the  skin  friction  is  determined  by 
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Quantities  Fc,  F^,  and  G  are  functions  of  the  Mach  number  at  the  edge  of  boundary -layer  and 
can  be  found  in  Ref.  67.  The  application  of  the  method  is  straightforward.  Some  typical  results 
are  shown  in  Fig.  13  along  with  other  theories68"70  and  experimental  data.71 


d.  Stratford -Beavers  Method 


Similar  to  the  foregoing  method,  semi-empirical  formulas  for  calculating  compressible 
turbulent  boundary-layers  were  developed  by  Stratford  and  Beavers.7-  The  formulas  for  the 
momentum  and  the  boundary-layer  displacement  thickness  are 
For  freestream  Reynolds  numbers  of  the  order  of  106, 
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For  freestream  Reynolds  numbers  of  the  order  of  10. 
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whore 


Res  = 


( 1  +  0.2M;)3‘w 


(64) 

(65 ) 

(66) 


1‘he  suffix  o  refers  to  stagnation.  The  quantity  w  is  the  exponent  in  the  viscosity-temperature 
relation  y  T00.  a  reasonable  value  for  cj  being  0.75,  except  at  very  high  temperatures,  where 

to  =  0.5. 

The  results  of  the  shape  factor  b*/0  found  by  Hqs.  (62)  and  (63)  compare  fairly  closely  to 
that  of  a  l/n^1  power  law.  where  n  =  7  and  c)  (see  Fig.  14). 


IV.  TRANSONIC  VISCOUS-INVISCID  INTERACTIONS 
1.  DESCRIPTION  OF  VISCOUS-INVISCID  TRANSONIC  FLOWS 

The  schematic  of  the  transonic  viscous-inviscid  interaction  over  an  airfoil  surface  is  repre¬ 
sented  in  Fig.  15.  The  embedded  supersonic  region  of  a  supercritical  flow  has  to  be  terminated 
by  a  shock  wave  to  bring  the  downstream  How  back  to  the  subsonic  state,  and  the  shock  foot  is 
smeared  into  a  series  of  compression  waves  as  a  result  of  viscous-inviseid  interactions.  The 
presence  of  an  adverse  pressure  gradient  due  to  the  compression  before  the  shock  usually  causes 
the  boundary  layer  to  separate  from  the  surface  for  either  laminar  or  turbulent  boundary  layers. 
However,  the  How  behavior  near  the  shock  differs  considerably,  depending  on  whether  the 
boundary  layer  is  laminar  or  turbulent  ahead  of  the  point  where  it  meets  the  shock. ^  The  pres¬ 
sure  rises  more  rapidly  for  turbulent  than  for  laminar  layers,^’4  while  the  displacement  thickness 
of  the  boundary  layer  increases  considerably  through  the  shock,  more  so  for  laminar  than  for 
turbulent  layers.-  Reynolds  number  has  a  strong  effect  on  the  interaction  in  the  case  of  laminar 
boundary  layers  but  almost  no  effect  for  turbulent  flows.-'4  In  any  event,  the  basic  features  of 
boundary-layer  thickening  and  pressure  rising  are  common  for  both  layers. 

Beyond  the  separation  point,  the  boundary-layer  How  is  divided  into  two  regions  by  the 
dividing  streamline.  The  How  above  this  streamline  includes  all  the  fluid  contained  in  the 
boundary  layer  just  upstream  of  separation.  Below  this  streamline,  the  How  yields  a  separated 
pattern  with  reversed  profile  near  the  wall.  Two  kinds  of  separated-flow  pattern  have  been  ob¬ 
served  by  Pearcey  et  al.^  for  turbulent  boundary  layers.  One  deals  with  a  local  separation  bubble 
caused  by  the  shock;  the  other  is  involved  with  rear  separation  and  depends  on  the  magnitude  of 
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tlie  pressure  gradient  approaching  the  trailing  edge  and  on  the  upstream  history  of  the  boundary 
layer. 

The  pressure  rises  continuously  in  the  subsonic  portion  although  its  gradient  decreases 
steadily  after  the  toe  of  the  shock. ^  This  feature  differs  from  that  of  a  supersonic  free  inter¬ 
action  in  which  the  incident  shock  is  reflected  as  an  expansion  fan  which  turns  the  flow  toward 
the  surface  and  thus  enables  the  reattachment  of  the  boundary  layer.  In  the  transonic  case, 
however,  the  shock  is  an  embedded  one  which  serves  as  part  of  the  compression  process  of  the 
flow  over  the  rear  of  the  airfoil;  whereas,  the  boundary  layer,  especially  the  laminar  one.  does 
not  have  sufficient  energy  for  reattachment  against  continuous  adverse  pressure  gradients  and, 
therefore,  remains  separated  all  the  way  toward  the  trailing  edge  (see  Fig.  16). 

The  mechanism  of  turbulent  viscous-inviscid  interaction  differs  from  that  of  laminar  inter¬ 
action.  In  the  turbulent  case,  when  the  flow  enters  the  strong  interaction  zone,  it  turns  away 
trom  the  surface  in  response  to  a  rapid  separation  bubble  growth,  triggered  from  the  toe  of  the 
shock; see  Ref.  5.  For  the  latter  case,  the  compression  before  the  shock  is  fairly  long  and  smooth 
and  involves  no  rapid  change  in  the  flow  angle  at  the  edge  of  the  boundary  layer;  see  Fig.  17. 

2.  STRONG  VERSUS  WEAK  INTERACTION  FORMULATIONS 

In  analysis  of  viscous  transonic  flows,  it  is  necessary  for  the  analytical  model  of  the  viscous 
system  to  have  the  capability  for  allowing  communication  of  positive  pressure  disturbance  from 
the  embedded  shock  wave  upstream  through  the  subsonic  portion  of  the  boundary  layer.  The 
surface  pressure  should  be  calculated  by  the  viscous  equations  rather  than  specified  by  the 
in  viscid  system.  The  communication  of  the  pressure  disturbance  through  the  viscous  layer  be¬ 
comes  possible  because  of  the  direct  influence  between  the  pressure  and  boundary  layer  thick¬ 
ness  in  a  subcritical  boundary  layer.  At  local  supersonic  speeds,  a  thickening  suberitieal 
boundary  layer  produces  a  rise  in  pressure  in  the  outer  stream  which  causes  further  thickening 
of  the  boundary  layer. ^ 

The  viscous  system  is  linked  to  the  inviscid  system  through  a  common  variable  by  which 
the  change  of  flow  properties  of  the  outer  flow  may  be  transmitted  to  the  inner  flow  or  vice 
versa.  It  reflects  the  physical  phenomenon  of  the  shock  wave'boundary-layer  interaction  process 
wherein  there  is  a  steep  pressure  rise  before  the  arrival  of  the  shock  wave  because  of  the  influence 
exerted  by  the  shock.  This  concept  of  the  new  viscous  system  has  been  explored  in  supersonic 
flows  by  Klineberg  and  Lees*'9  and  Lees  and  Reeves^  and  in  transonic  flows  by  Klineberg  and 
Steger,^  and  Tai. 'OJS 

The  new  system  model  is  referred  to  as  the  strong  interaction  formulation.  Other  formula¬ 
tions,  such  as  the  usual  boundary-layer  system  by  which  the  viscous  effect  is  accounted  for  by 
the  displacement  correction  or  analytical  techniques  which  do  not  involve  flow  separation,  are 
referred  to  as  the  weak  interaction.  The  strong  interaction  formulation  may  be  applied  to 
attached,  as  well  as  separated,  boundary  layers.  When  applied  to  attached  flows,  the  boundary 
layer  eventually  separates.  If  the  same  flow  were  to  be  treated  by  the  weak  interaction  system. 
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numerical  experiments  indicate  that  the  boundary  layer  would  remain  attached  until  shock  jump 
was  encountered.  This  gives  further  indication  that  the  strong  interaction  system  is  more  suitable 
in  simulating  the  shock  wave  boundary-layer  interaction  process  than  are  the  usual  boundary- 
layer  approaches.  On  the  other  hand,  the  use  of  weak  interaction  formulation  is  preferred  in  the 
forward  portion  of  the  airfoil  where  the  viscous-inviscid  interaction  is  presumably  weak. 

3.  CALCULATION  OF  WEAK  INTERACTIONS 

a.  Boundary  Layer  Correction 

As  in  the  usual  procedure,  the  boundary-layer  quantities  are  calculated  based  on  the  speci¬ 
fied  in  viscid  pressure  distribution,  and  the  inviseid  solution  is  updated  based  on  the  surface  aug¬ 
mented  by  the  boundary-layer  displacement  thickness.*  In  so  doing,  the  two  Hows  are  con¬ 
nected  at  the  locus  of  the  boundary-layer  displacement  thickness.  The  latter  serves  as  a  '‘stream¬ 
line"  through  which  mass  transfer  or  more  importantly,  direct  communication  between  the  outer 
inviseid  and  inner  boundary-layer  Hows  is  prohibited.  The  inviseid  solution  is  then  repeated 
until  a  converged  solution  is  obtained  upon  satisfying  converging  criterion  after  several  iterations. 

Existing  transonic  codes  by  Bauer  et  al.,74  and  Carlson7^  determine  the  viscous  effects  by 
such  a  procedure.  The  Nash-Macdonald  method67  for  turbulent  boundary  layers  is  used  in  both 
codes  (weak  interaction  in  the  laminar  case  is  simply  neglected).  In  the  case  of  transonic  cascade 
flows,  Cdiebe  coupled  an  inviseid  flow  solution  by  a  finite-difference  method  and  a  viscous  flow 
solution  by  the  Stratford-Beavers  method.77  Fig.  18  shows  typical  results  of  the  boundary- 
layer  effect  on  the  blade  of  a  cascade  flow  evaluated  by  Gliebe.18 

b.  Analytical  Method 

Mason  and  Inger14  have  developed  an  analytical  method  for  calculating  the  transonic 
viscous-inviscid  interaction  involving  a  weak  normal  shock  wave  impinging  on  unseparated  turbu¬ 
lent  boundary  layers.  The  interaction  field  is  considered  as  a  multisubregion  problem  and  the 
resulting  viscous  transonic  tlow  system  is  solved  by  the  Fourier  transformation  method.  The 
theoretical  model  is  valid  for  weak  shock  conditions  with  the  Mach  number  ahead  of  the  shock 
not  very  near  unity.  The  schematic  of  the  model  which  is  shown  in  Fig.  19  includes  outer  mixed 
tlow  and  inner  rotational  disturbance  flow  with  a  Lighthill  sublayer. 

c.  Asymptotic  Expansion  Method 

Melnik  and  Grossman, ^  on  the  other  hand,  formulated  an  asymptotic  expansion  method 
for  analyzing  the  interaction  of  a  normal  shock  wave  with  an  unseparated  turbulent  boundary 

♦Another  way  of  accounting  for  the  boundary-layer  effect  is  to  determine  an  equivalent  source  distribution 
o  =  d(Utt6  *)/ds  as  suggested  by  Lighthill.76 


layer  on  a  Hat  surface  at  transonic  speeds.  The  interaction  involves  an  outer,  inviscid  rotational 
layer,  a  constant  shear  stress  layer  near  the  wall,  and  a  blended  layer  between  the  two.  Fig.  20 
shows  typical  results  of  the  asymptotic  method  compared  with  those  of  the  Manson-Inger 
scheme.14 


4.  CALCULATION  OF  STRONG  INTERACTIONS 
a.  Laminar  Viscous-lnviscid  Interaction 


The  calculation  is  illustrated  by  using  the  integral  boundary  layer  method  coupled  with  an 
inviscid  solution  by  the  method  of  integral  relations.10  The  coordinate  system  for  the  coupling 
is  shown  in  Fig.  2  1 .  The  ordinary  differential  equations  for  the  inviscid  external  How,  reduced 
by  the  method  of  integral  relations,  are  (for  simplicity  only  those  at  the  edge  of  the  boundary 
layer  are  written  here) 
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and  those  for  the  viscous  system10  are 
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The  viscous  system  is  coupled  directly  to  the  inviscid  system  by  the  induced  angle  of  invis¬ 
cid  streamline  at  the  edge  of  the  boundary  layer 


0  =  sin" 1 


Ve  (invis.) 
Meac  (viscous) 


<(V7) 


where  Ve  is  calculated  by  the  inviscid  system,  while  the  velocity  magnitude  Mcae  is  obtained  by 
the  viscous  system.  The  variable  0  is  a  common  variable  for  both  inviscid  and  viscous  systems 
and,  therefore,  governs  the  viscous-inviscid  interaction  process.  Since  the  two  Hows  are  coupled 
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by  the  inviscid  streamline  angle  rather  than  the  streamline  itself,  inass  transfer  between  the  outer 
in  viscid  and  inner  boundary-layer  Hows  is  allowed  in  accordance  with  the  continuity  equations 
by  which  the  variable  0  is  introduced  in  the  viscous  system 


d6* 

ds 


tan  0  +  (6  -  6*) 


(68) 


The  value  of  0.  which  is  determined  by  the  viscous-inviscid  interaction,  has  a  direct  bearing  on 
the  growth  of  the  boundary  layer. 

The  strong  interaction  system  is  applied  at  some  distance  preceding  the  inviscid  shock  loca¬ 
tion.  The  exact  initial  location  of  the  strong  interaction  region  is  determined  by  an  iterative 
process. 10 

Results  of  calculations  at  supercritical  freestream  Mach  numbers  are  presented  for  both  a 
6  percent  circular  arc  and  an  NACA  0015  airfoil. 

Fig,  22  gives  the  boundary-layer  displacement  thickness  throughout  a  6  percent  circular 
arc  airfoil  at  M^  =  0.868  and  Re^  =  6.9  x  104.  The  thickening  of  the  boundary  layer  in  the 
forward  portion  follows  a  similar  trend  as  that  found  by  Schubauer,  using  the  Karman-Polhausen 
method77 ;  however,  the  strong  interaction  calculation  gives  a  far  more  realistic  6  distribution 
pattern  in  the  rear  portion.  The  boundary  layer  is  practically  of  the  Blasius  type  (a  =  1 .857 )  in 
the  leading-edge  region  and  varies  slightly  throughout  the  forward  portion  of  the  airfoil.  It  re¬ 
mains  unseparated  through  the  embedded  supersonic  region  although  the  viscous-inviscid  inter¬ 
action  has  been  strong  since  x  =  0.37.  The  separation  point  is  found  when  a  =  0  which  corres¬ 
ponds  to  zero  shearing  stress  at  the  wall. 

The  boundary  layer  keeps  separated  over  the  rear  of  the  airfoil  where  small  adverse  pressure 
gradients  are  generated  by  continuous  compression  of  the  outer  subsonic  flow.  This  is  a  physical 
feature  of  the  transonic  viscous-inviscid  interaction  since  by  compression,  the  How  ought  to 
return  almost  to  the  freestream  value  downstream.  After  leaving  the  trailing  edge,  the  viscous 
wake  contains  a  reversed  How  for  a  distance  and  then  turns  to  forward  tlow  downstream.  The 
location  of  the  rear  stagnation  point  agrees  well  with  that  found  bv  Klineberg  and  Steger^  under 
similar  How  conditions. 

Theoretical  results  of  pressure  distribution  for  the  6  percent  circular  arc  airfoil  compare 
very  well  with  recent  laminar  experimental  data  of  Collins  and  Krupp78  as  presented  in  Fig.  23. 
not  only  in  the  pressure  distribution  but  also  in  the  separation  point.  The  small  difference  in 
freestream  Mach  number  between  theory  and  experiment  was  deliberately  selected  to  offset 
wind-tunnel  interference  effects.79  Fig.  23  also  shows  the  inviscid  isentropic  and  nonisentropic 
solutions  obtained  by  using  the  method  of  integral  relations  with  three  different  shock  condi¬ 
tions.  Note  that  the  shock  foot  is  smeared  as  a  result  of  the  strong  viscous-inviscid  interaction. 
The  inviscid  result  with  AS  >  0  and  0  =  75  degrees  simulates  the  How  best  in  the  region  near 
the  shock ;  however,  it  still  deviates  from  the  viscous  data  in  the  trailing  edge.  The  difference  is 
obviously  attributable  to  the  thickening  of  the  laminar  boundary  layer  toward  the  trailing  edge. 
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On  the  other  hand,  comparisons  between  theory  and  experiment  for  the  NASA  0015  air¬ 
foil  were  jeopardized  because  of  lack  of  available  laminar  test  data.  Measurements  available  for 
this  case  were  taken  from  Graham  et  al.H0  at  Reynolds  numbers  on  the  order  of  2  x  I06;  the 
flows  were,  therefore,  turbulent.  The  calculated  laminar  result  is  in  good  agreement  with  the 
experimental  data  in  the  forward  portion  of  the  airfoil;  however,  there  is  considerable  dis¬ 
crepancy  in  pressure  near  the  trailing  edge;  see  big.  24.  The  deviation  may  be  attributed  to  the 
extended  separation  of  the  laminar  boundary  layer,  which  would  be  unlikely  to  occur  if  the 

boundary  layer  were  turbulent.  Since  the  method  of  using  reversed-flow  velocity  profile  was 

x  1 

originally  attempted  for  application  to  cases  with  small  adverse  pressure  gradients01 ,  this  dis¬ 
crepancy  also  raised  questions  regarding  the  validity  of  its  use  in  separated  boundary  layers  with 
strong  adverse  pressure  gradients. 

Finally,  Fig.  25  shows  the  overall  How  Held  of  the  NACA  0015  airfoil  at  =  0.729, 
a  =  4  degrees,  and  Re^  =  0.145  x  106,  based  on  the  chord  length.  The  calculated  boundary 
layer  is  fairly  thin  because  of  the  high  Reynolds  number  used.  The  separated  boundary  layer  on 
the  lower  side  is  thicker  than  that  on  the  upper  side  because  the  latter  is  developed  from  a  ver> 
thin  attached  layer  in  an  accelerating  external  flow  prior  to  the  shock  wave.  In  any  case,  the 
basic  features  of  the  transonic  viscous-in  viscid  interaction  are  clearly  seen.  The  compression 
waves  in  the  region  near  the  shock  are  represented  by  the  Mach  lines.  Being  normal  to  the  edge 
of  the  boundary  layer,  the  Mach  one  line  represents  the  end  of  a  smearing  shock  toe  rather  than 
the  whole  shock.  The  mechanism  of  the  compression  has  been  designated  in  the  literature  as 
the  “lambda  shock,"  which  serves  as  part  of  the  decelerating  process  of  the  How  over  the  rear  of 
the  airfoil,  where  the  laminar  boundary  layer  does  not  have  sufficient  energy  for  reattachment 
against  continuous  adverse  pressure  gradients  and,  therefore,  remains  separated  all  the  way  to 
the  trading  edge. 

It  is  also  interesting  to  note  that  the  viscous-inviscid  interaction  process  terminates  the  em¬ 
bedded  supersonic  How  earlier  than  does  the  inviscid  shock.  The  result  is  not  surprising  since  the 
interaction  process  has  weakened  the  strength  of  the  shock  wave;  in  turn,  it  moves  the  shock 
forward.  However,  because  of  the  change  of  the  local  velocity  distribution,  the  height  of 
the  supersonic  pocket  determined  by  the  theory  is  slightly  higher  than  that  obtained  by  the 
inviscid  solution;  it  remains  to  be  verified  by  appropriate  experiment  and  other  valid  analysis. 

The  integral  boundary  layer  method  can  also  be  coupled  with  a  finite  difference  scheme 
carried  out  by  Klineberg  and  Steger  8  There  the  V  values  were  prescribed  rather  than  calcu¬ 
lated.  Although  the  prescribed  Ve  values  were  improved  successively,  the  viscous  and  inviscid 
flows  had  to  be  treated  separately  even  for  the  strong  interaction.  More  important ly.  in  the 
indirect  coupling  as  such,  it  would  be  difficult  to  allow  direct  communication,  including  the  mass 
transfer,  between  the  outer  inviscid  flow  and  the  boundary  layer  flow. 
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b.  Turbulent  Viscous-lnviscid  Interaction 


Because  of  the  convenience  for  direct  coupling  of  the  outer  inviscid  and  viscous  flows,  the 
method  of  integral  relation  is  used  again  for  the  inviscid  solution  in  the  illustration  for  calculating 
the  turbulent  viscous-inviscid  interaction.1**5  The  coordinate  system  and  inviscid  tlow  equations 
are  the  same  as  for  the  laminar  case;  Eqs.  (30)  and  (31 ).  The  equations  for  the  viscous  system 
are  derived  by  adding  the  continuity  equation  to  the  original  Kuhn-Nielsen  method,12  i.e.. 


a  rt  -  d«  _  _ 
Hi  J lldT? '  l,eHF_ ' v< 

o 

I  he  resulting  viscous  system  islri 


(69) 


(70) 


where  Hy  and  Oj  are  presented  in  Ref.  15.  The  system  is  valid  for  both  attached  and  separated 
Hows. 

Similar  to  the  laminar  analysis,10  the  viscous  system  is  coupled  directly  to  the  inviscid  sys¬ 
tem  by  the  induced  angle  of  the  inviscid  streamline  at  the  edge  of  the  boundary  layer,  which 
allows  mass  transfer  between  the  two  flows  in  accordance  with  the  continuity  equation.  Eq.  (69). 
The  strong  interaction  system  is  applied  at  some  distance  before  the  inviscid  shock  wave  location. 
As  opposed  to  the  laminar  case,  the  tlow  angle  at  the  edge  of  the  boundary  layer  is  adjusted. 

The  reason  for  the  adjustment  is  that  in  the  turbulent  case,  when  the  tlow  enters  the  strong 
interaction  zone,  it  turns  away  from  the  surface  in  response  to  a  rapid  separation  bubble  growth, 
triggered  from  the  toe  of  the  shock ;  see  Ref.  5.  The  boundary  layer  is  attached  at  the  start  of 
the  strong  interaction  region;  however,  it  usually  separates  in  a  short  distance.  Computation  is 
then  switched  to  the  subroutines  based  on  separated  velocity  profiles  of  the  strong  interaction 
system. 

Results  of  calculations  at  supercritical  freest  ream  Mach  numbers  are  presented  for  a  10- 
percent  thickness  bump  and  an  NACA  0015  airfoil  at  an  angle  of  attack.  The  bump  is  basically 
a  10-percent  circular  arc,  faired  with  cosine  curves  at  both  ends.  big.  2(>  gives  the  typical  velocity 
profile  at  different  stations  over  the  10-percent  thickness  bump  at  M^  =  0.7325  and  Rc^  -  1 .75 
x  I06.  The  profile  at  x  =  0.2  is  the  initial  profile,  calculated  by  using  Schlichting's  skin  friction 


formula  based  on  a  power  law  distribution.^  (iood  comparison  between  the  calculated  and 
available  measured  profiles  of  Alber  et  al.<>(>  at  x  =  0.5  and  0.875  indicates  the  adequacy  ol  the 
theoretical  approach. 

fig.  27  shows  the  distribution  of  the  friction  velocity  over  the  same  bump.  Note  that  the 
friction  velocity  in  the  strong  interaction  region  is  calculated  in  the  analysis,  rather  than  pre¬ 
scribed  as  in  Ref.  12.  The  How  separates  in  a  fairly  short  distance  after  it  enters  the  strong  inter¬ 
action  zone  influenced  by  the  shock.  The  friction  velocity  remains  almost  constant  downstream 
of  the  separation  point  and  then  increases  gradually  toward  reattachinent .  The  calculated  pres¬ 
sure  distribution  compares  fairly  well  with  recent  turbulent  experimental  data  ot  Alber  et  al  06  as 
presented  in  Fig.  28.  Also  plotted  are  the  inviscid  solutions  obtained  by  using  the  method  ol 
integral  relations,  with  and  without  entropy  change  across  the  shock,  file  agreement  between 
theoretical  and  experimental  results  is  good  both  for  pressure  distribution  and  separation  point. 
Note  that  the  shock  foot  is  smeared  as  a  result  of  the  strong  viscous-in  viscid  interaction.  1  he 
compression  starts  at  x  =  0.575  where  the  strong  interaction  begins.  The  turbulent  boundary 
layer  is  attached  throughout  the  supersonic  region  and  the  strong  interaction  starts  in  front  ot 
the  shock  wave.  Flow  separation  takes  place  downstream  of  the  shock  wave  at  a  peak  Mach 
number  of  1.14.  The  trend  is  consistent  with  that  found  experimentally  by  Alber  et  al.°°  lor 
cases  with  Mp  <  1 .32.  The  turbulent  boundary  layer  reattaches  downstream  of  the  bump. 

The  difference  between  the  theoretical  and  experimental  pressures  in  the  rear  of  the  bump  is 
attributed  to  the  insufficient  damping  effect  in  the  inner  layer  eddy  viscosity  model.  There  is 
also  reason  to  believe  that  it  may  be  due  to  overprediction  of  the  flow  by  the  strong  interaction 
equations. 

The  flow  over  an  NACA  0015  airfoil  at  =  0.729  and  ot  -  4  degrees  has  been  investi¬ 
gated  as  an  example  for  the  lifting  case.  This  particular  flow  condition  has  been  calculated  in 
the  laminar  analysis.10  There  the  agreement  between  the  laminar  theory  and  the  experiment 
was  good  except  in  the  rear  of  the  airfoil,  where  the  flow  could  have  become  turbulent .  There¬ 
fore,  the  turbulent  viscous-inviscid  interaction  model  is  employed  for  flows  in  the  rear  portion 
of  the  airfoil.  The  transition  zone  is  assumed  to  he  short  enough  so  that  the  continuity  of 
boundary  layer  properties  may  be  maintained.  The  corresponding  initial  turbulent  shear  velocity 
is  then  estimated  based  on  similar  velocity  profiles  for  both  laminar  and  turbulent  boundary 
layers.  The  weak  interaction  system  is  used  in  actual  computation  after  the  strong  interaction 
formulation  fails  to  produce  a  converged  solution.  Calculated  results  are  given  in  Fig.  29  along 
with  the  previous  laminar  solution10  and  the  experimental  data  of  Graham  et  al*° 

As  indicated  in  Fig.  29,  the  pressure  values  over  the  rear  region  of  the  airfoil  calculated  by 
the  present  turbulent  theory  compare  fairly  well  with  the  experimental  data.  Ref.  80.  which  were 
measured  at  Reynolds  numbers  on  the  order  of  2  x  I0(\  The  flow  over  the  upper  surface  lias 
undergone  a  strong  viscous-inviscid  interaction  in  the  vicinity  of  the  shock  wave  between  \  = 

0.34  and  0.6,  but  it  seems  to  have  returned  to  a  weak  interaction  process  after  x  =  0.6;  even  the 
boundary  layer  is  still  separated  thereafter.  The  boundary  layer  remains  separated  all  the  way  to 
the  trailing  edge.  It  follows  a  similar  flow  pattern  to  that  observed  by  Pcarcey  et  al.> 
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5.  VISCOUS-INVISCID  INTERACTION  IN  SEPARATED  FLOW  REGION 


Although  tlk'  strong  mk'Mclioii  lormulation  is  use*  I  ill  ami  necessary  lor  Ilk*  shock  wave 
region,  it  sutlers  a  problem  of  accuracy  ami  numerical  stability  in  the  separated  flow  region  I  he 
reason  for  this  is  believed  to  be  attributed  10  the  basic  feature  ol  the  strong  interaction  system 
that  the  pressure  is  calculated  simultaneously  with  the  boundary  layer  thickness  I  he  results  of 
the  integral  methods  are  therefore,  much  more  sensitive  to  the  velocity  profiles  used  tor  either 
laminar  or  turbulent  boundary  layers.  Any  inadequacies  in  velocity  protiles  employed  will  be 
magnified  in  the  strong  interaction  system,  as  opposed  to  being  suppressed  in  the  usual  boundary 
layer  calculations. 

In  general,  the  velocity  profiles  for  separated  Hows  are  more  difficult  to  define  than  for 
attached  boundary  layers.  In  the  foregoing  laminar  case,  the  reversed-flow  profile  was  originally 
attempted  for  application  to  cases  with  small  adverse  pressure  gradients  s|  Its  use  in  separated 
boundary  layers  with  strong  adverse  pressure  gradients  introduced  some  discrepancy  m  agree¬ 
ment  between  the  theory  and  the  experiment  in  the  trailing  edge  region.10 

In  the  case  of  turbulent  Hows,  on  the  other  hand,  a  numerical  difficulty  is  encountered  in 
using  the  strong  interaction  system.  Very  fine  adjustment  of  friction  velocity  u  lias  to  be  made 

for  achieving  the  convergence  of  the  flow  downstream  of  the  separation  point.1  N  see  Fig.  30 
Because  of  the  convergence  problem,  it  is  proper  to  use  t lie  weak  interaction  system  (pressure 
prescribed)  rather  than  the  strong  interaction  formulation  near  and  aft  the  separation  point. 

These  arguments  suggest  that  in  the  separated  flow  region,  a  weak  interaction  system  might 
provide  more  reliable  results  than  the  strong  interaction  formulation.  I  he  point  of  switching 
from  the  strong  to  the  weak  interaction  equations  has  to  be  determined  by  numerical 
experiments. 


V.  APPLICATION  TO  CASCADE  FLOWS 

The  flow  in  a  cascade  at  transonic  speeds  is  far  more  complex  than  the  steady  .  two-dimen¬ 
sional  flow  past  an  airfoil  in  free  air.  A  comprehensive  review  of  the  tlow  m  transonic  com¬ 
pressors  has  been  given  by  Kerrebrock  ***  To  apply  the  previous  techniques  to  cascade  Hows 
involving  transonic  viscous-inviscid  interactions,  assumptions  that  are  necessary  are  that  the  tlow 
m  a  cascade  has  axisymmetric  stream  surfaces  and  that  the  flow  on  these  surfaces  can  be  treated 
two  dimensionally.  Such  Hows  very  nearly  exist  in  cascades  where  the  blade  height  is  small 
compared  to  the  radius  from  the  centerline  of  the  cascade.  I  hey  are  called  quasi  two-dimen¬ 
sional  flows. I  he  schematic  of  a  quasi  two-dimensional  flow  m  an  axisymmetric  stream  sur¬ 
face  is  depicted  in  Fig.  3  1.  I  he  approximation  represents  an  ideali/ed  situation  lot  the  complex 
unsteady  flow  in  transonic  compressors.  It  is.  nevertheless,  useful  foi  qualitative  investigation  ol 
the  influence  of  various  design  parameters. 
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The  governing  equation  for  a  quasi  two-dimensional  compressible  in  viscid  How  in  an 

•  xs 

axisyniinctric  stream  surface  represented  by  coordinates  (£,  0  )  is  • 
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where 


ca  =  constant  rotation  of  the  cascade 
R  =  distance  measured  from  the  axis 
h  =  blade  height 


The  stream  function  \p  is  defined  by 


=  phV 


if  =  -  PhRU 


The  boundary  conditions  are  q  =  0  on  the  blade  surfaces. 

Eq.  (71 )  can  be  solved  by  the  finite  difference-relaxation  method,  similar  to  the  procedure 
employed  in  a  two-dimensional  potential  flow  about  an  airfoil. 

To  facilitate  calculating  the  strong  viscous-in  viscid  interactions  and  to  account  for  the  non- 
isen tropic  effect  of  the  embedded  shock,  it  was  suggested  in  the  foregoing  discussions  that  the 
hybrid  method  is  mostly  suitable  for  the  purpose.  That  is.  the  How  in  the  shock  region  is  solved 
by  the  method  of  integral  relations  with  the  initial  and  boundary  conditions  provided  in  the 
finite  difference-relaxation  approach.  To  implement  the  method  of  integral  relations,  the 
governing  equations  for  the  quasi  two-dimensional  channel  How  between  the  cascade  blade 
surfaces  are  written  in  primitive  form  as 
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With  u  boundary  condition  of  zero  normal  velocity  on  the  surface,  Lqs.  (73)  through  (76) 
can  be  solved  by  the  method  of  integral  relations  with  the  aid  of  the  N-2-strip  scheme  4“  The 
final  solution  is  achieved  by  updating  solutions  between  the  finite  difference  method  tor  the 
overall  flow  and  the  method  of  integral  relations  for  the  shock  region. 

The  in  viscid  solution  so  obtained  is  convenient  to  be  coupled  with  either  laminar  or  turbu¬ 
lent  viscous  systems  for  strong  interaction  between  the  shock  wave  and  the  boundary  layer. 

A  weak  interaction  formulation  should  be  applied  for  conventional  boundary  layer  regions, 
including  a  separated  flow  /one  downstream  of  the  shock  wave. 
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(Illustration  is  based  on  Ref.  32) 
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Figure  4  -  Extension  of  Freestream  Boundary  to  ‘infinity”  for  Flowfield  Solved  by 

the  Method  of  Integral  Relations 
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F  igure  5  Pressure  Distribution  on  an  N  AC  A  64A410  Airfoil  at  =  0.72  and  q  =  4 


FigureH  -  Pressure  Distribution  Over  an  NACA  64A006  Airfoil  at  a  =  0 
degree  and  =  O.S75  CaJculated  by  the  Finite  Element  Method 


Figure  9  Comparison  Between  the  Numerical  Results  Obtained  by 
the  Finite  Element  Method  and  the  Finite  Difference  Method 
for  Transonic  Small  Disturbance  Equation 
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Figure  10  —  Theoretical  Pressure  Distributions  Around  an  NACA 
64A410  Airfoil  at  -  0.72,  and  a  =  0  Degrees 
(Illustration  based  on  Ref.  58) 


Figure  1 1  Results  of  Adiabatic  Laminar  Boundary -Layer  Shock -Wave 
Interaction  by  the  Klineberg-Lees  Method  Over  Flat  Plate 
at  Supersonic  Speed 
(Illustration  based  on  Ref.  59) 
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Figure  12  Results  of  an  Adiabatic  Turbulent  Boundary  Layer  Flow  Past  a  Bump 
at  Transonic  Speed  by  Kuhn-Nielsen  Method 
(Illustration  based  on  Ref  12) 


Figure  M  Momentum  Thickness  of  a  Compressible  Turbulent  Boundan  La\cr  Calculated 
by  Nash-Macdonald  Method  Together  with  Others 
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Figure  14  Profile  Shape  Parameter  h*;U  of  a  Compressible  Turbulent  Boundary  Layer 
Calculated  by  Stratlord-Beavers  Method 
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Figure  15  Transonic  Viscoiis-lnviscid  Interaction  Overall  Airfoil 
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Figure  16  -  Schematic  of  Supersonic  and  Transonic  Viscous-In  viscid  Interactions 
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Figure  17  Comparison  of  Laminar  and  Turbulent  Viscous-In  viscid  Interactions 
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Figure  21  -  Coordinate  Systems  and  Flow  Regions  for  Calculation  of 
Viscous-Inviscid  Transonic  Flow 


Figure  22  Boundary-Layer  Displacement  Thickness  for  a  6rr  Circular  Arc 
at  Mo,,  =  0.868  and  Re^  =  6.9  x  I04 
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Figure  23  -  Pressure  Distribution  Over  a  6 % 
Circular  Arc  at  M ^  -  0.868  and  a  =  0° 
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Figure  24  Pressure  Distribution  on  an  NACA  0015  Airfoil  at  =  0.729  and  a  =  4° 
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(NORMALIZED  BY  V, 


Figure  29  -  Pressure  Distribution  on  an  NACA  0015  Airfoil  at  “  0.729 

and  a  =  4  Degrees 
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Figure  50  Friction  Velocity  Gradients  Downstream  of  the  Separation  Point 

<  Ref.  15) 
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RADIAL  DIRECTION 


Figure  31  -  Quasi  Two-Dimensional  Flow  on  Axisymmetric  Stream  Surface 
(Illustration  based  on  Ref.  85) 
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